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Abstract. 

Numerical simulations show that microscopic rod-like bodies suspended in a 
turbulent flow tend to align with the vorticity vector, rather than with the dominant 



' eignevector of the strain-rate tensor. This paper investigates an analytically solvable 

^ . limit of a model for alignment in a random velocity field with isotropic statistics. The 

I vorticity varies very slowly and the isotropic random flow is equivalent to a pure strain 

^ ' with statistics which are axisymmetric about the direction of the vorticity. We analyse 

^ , the alignment in a weakly fluctuating uniaxial strain field, as a function of the product 

' ' of the strain relaxation time Tg and the angular velocity ui about the vorticity axis. We 

^~>i find that when luTs S> 1, the rods are predominantly either perpendicular or parallel 

'" ^ to the vorticity. 

> 

00 ■ 1. Introduction 

\o 

\^ [ Microscopic rod-like bodies suspended in a fluid flow rotate in response to the velocity 

l/^ I gradient of the flow. This introduces a degree of order in the orientation of a suspension 



of particles which can influence its optical or rheological properties. The equation of 
motion for the orientation of microscopic ellipsoidal particles was obtained by Jeffery 
[1] . The implications of this equation of motion for the orientation have been considered 
^ I by numerous authors: for example [2] discusses the motion of general axisymmetric 
particles, [3] considers the role of Brownian motion, |H [5] discuss the alignment 
fields in (respectively) regular and chaotic flows, and [6l [7] are recent experimental 
contributions which give an extensive list of references. There are, however, still 
aspects which are not thoroughly understood. One surprising observation (based 
upon direct numerical simulation (DNS) studies of Navier-Stokes turbulence) is that in 
isotropic fully-developed turbulence, rod-like particles show significant alignment with 
the vorticity vector, but negligible alignment with the principal strain axis [8]. This 
was given a qualitative explanation in [8], but it is desirable to have a model for this 
surprising effect which can be analysed quantitatively. 

This paper considers an exactly solvable model for the alignment of rods with 
vorticity. The formulation of this model was motivated by observations about the 
velocity gradient field of turbulence. It has been observed that the fluctuations of the 
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vorticity vector decay much more slowly than fluctuations of the rate of strain: [S] shows 
evidence that the correlation functions of strain and vorticity both show approximately 
exponential decay, with decay times Tg ~ 2.3tk and ~ 7.2rK respectively, where tk 
is the Kolmogorov timescale of the turbulence. Similar results were reported earlier by 
Girimaji and Pope [9] and Brunk, Koch and Lion [10]. This observation suggests that 
it may be helpful to consider the limit as — )■ oo, that is the limit where the vorticity 
is frozen, in order to explain the observed alignment. 

We use an Ornstein-Uhlenbeck process to model fluctuations of the velocity 
gradient, and consider the limit where the vorticity evolves very slowly. This model 
is solved exactly in the limit where the strain which occurs over the timescale Ts is 
small. The alignment of the rod direction n and the direction of the vorticity vector 
can be described by computing the probability density function (PDF) of 2; = n • e^^^. 
We find that in these limits the PDF of z, denoted by P{z), can be computed exactly. 
This analytically solvable model has a single dimensionless parameter, ( = ut^, where 
u is the angular velocity of rotation about the vorticity vector. We find that when 
C ^ 1, the probability density has two sharp peaks, one at 2; = ±1 (indicating perfect 
alignment with vorticity), the other at 2 = (implying that the rods are perpendicular 
to the vorticity). In the limit as ( ^ 00, the peak at z = ±1 is higher than at z = 0, but 
it is also narrower, with both peaks containing a finite probability. (Throughout this 
paper, {X) is the expectation value of X, and we use P{X) to denote its probability 
density function). 

Section [2] discusses the model which will be solved: the equations of motion for a 
microscopic rod are considered in section 12.11 and the Ornstein-Uhlenbeck model for 
the velocity gradient of an isotropic random flow is described in section 12.21 Section 
[3] discusses a transformation of the equation of motion in which the isotropic velocity 
gradient is replaced by a pure strain field which is axisymmetric about the direction of 
the vorticity vector, and it discusses the parametrisation of such axisymmetric random 
strain fields. Section H] considers the general solution for alignment of rod-like particles 
in axisymmetric strain fields, before specialising to the solution of the model developed 
in section [3l Section [5] summarises our conclusions. The analysis in section H] is closely 
related to recent work by Vincenzi pi], who analysed the alignment of ellipsoidal 
particles in an axisymmetric Kraichnan-Batchelor model. 

2. Equations of motion 

2.1. Non-linear and linear equations of motion for rods 

We consider microscopic objects advected in a fluid with velocity field u{r,t). The 
objects are assumed to be neutrally buoyant, and smaller than any lengthscale 
characterising the fluid, but sufficiently large that their Brownian motion need not 
be considered. The motion of the body is described by the position of its centre, r{t), 
and the direction of a unit vector aligned with its axis, n(t). The centre of the body is 
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assumed to be advected by the fluid flow: r = u{r, t). The motion of the unit vector n 
defining the axis of symmetry is determined by elements of the velocity gradient tensor, 
evaluated at the centre of the body: 

Mt) = ^^ir{t)^t) (1) 

where r{t) is the advected particle trajectory. The equation of motion of the director 
vector of a microscopic rod-like body is [I] 

^ = A(t)n-(n.A(t)n)n. (2) 

We assume the flow is incompressible, so that Z]f=i = 0. This tensor can 
be decomposed into a symmetric part S, which is termed the strain rate, and an 
antisymmetric part f2, which is the vorticity tensor: 

A = S + f2 , = s , = -fl . (3) 

If the velocity gradient matrix were constant in time, the equation of motion ([2]) 
would imply that the vector n would become aligned with the eigenvector corresponding 
to its largest eigenvalue. However, numerical simulations of equation ([2]) for velocity 
flelds of fully developed turbulence show a different, and unexpected, phenomenon [8]. 
It is found that the direction vector n has negligible correlation with the dominant strain 
eigenvector, but that it does have a quite pronounced correlation with the vorticity 
vector, uj. 

Our analysis of the alignment due to the motion ([2]) will use an observation due to 
Szeri [12]: the non-linear equation ([2]) can be solved by considering a companion linear 
equation for a vector cc(t), which evolves under the action of a monodromy matrix M(t): 

x{t) = M(t)a;(0) , ^ = A(t) M (4) 

where the initial conditions are M(0) = I (the identity matrix) and x{0) = n(0). The 
solution to (E]) is obtained by normahsing the solution of @: 

- rI • 

The advantage of this approach is that it is easier to solve the linear equation (jlj) than 
the non- linear equation ([2]). 

2.2. Ornstein- Uhlenbeck model for velocity gradients in isotropic flows 

In this section we describe a simple stochastic model for the matrix A(t) in isotropic 
random flows. A version of this model was used by Vincenzi et al [13] , and its structure 
is suggested by the observations in [10]. The model was also considered in [8], which 
gave a detailed account of its implementation. Here we give a brief summary. 

It is known that the elements of S and Q fluctuate randomly about zero, 
with different timescales and respectively. Their correlation functions are well 
approximated by exponential functions. This suggests modelling the elements of S and 
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O by Ornstein-Uhlenbeck processes The three independent components of the 

vorticity will be modelled by: 

1 



Q = Q,^ + J2D,r]ij{t) (6) 

T 
' V 

where the riij{t) are independent white- noise signals, satisfying 

m) = o, mv{t')) = 6{t-t'). (7) 

This model predicts that the correlation function of Qij is exponential [HI [T2] : 

(a,(ti)f],,(t2)) = Dvrvexp(-|ti - tal/T-v) • (8) 

The components of the strain-rate matrix are generated by a further six Ornstein- 
Uhlenbeck processes, with a different correlation time Ts. The off-diagonal elements 
are generated by a process of the same form as (Q, with the diffusion coefficient 
in ([7]) replaced by Dg. The diagonal elements of the strain-rate matrix must satisfy 
J2i=i Sii = 0, which is the incompressibility condition, V ■ u = 0. This constraint is 
satisfied by the solution of the following Ornstein-Uhlenbeck equations 

3 



Sii — Sii + \/2Df] 



(9) 



The elements flij and Sij generated by these processes are statistically independent, 
apart from the constraint that J2i Su = 0. The variances of the off-diagonal, diagonal 
and vorticity elements are respectively denoted (S^), (S'J) and (fi^), and are related to 
the relaxation times and diffusion rates by {S"^) = DsTs, {S^) = ^D^Tg and {Q"^) = D^t^. 
The requirement that the statistics of the model are invariant under rotations (so that 
it describes a velocity gradient with isotropic statistics) gives = 2Ds, so that this 
model has four parameters: Tg, t^, Dg and D^. Note that the diffusion coefficients have 
dimension [D] = T~^, implying that the model has three independent dimensionless 
parameters. In the following we consider the limit as t^/ts — oo, so that the vorticity 
is frozen, with angular velocity u. We also assume that DsT^ ^ 1, so that the strain 
fluctuations are small. This leaves one dimensionless parameter, which we will take to 
be C = i^Ts- 



3. Transformation to an axisymmetric pure strain model 

3.1. The frozen vorticity limit 

In this section we consider the alignment of rod-like particles in an isotropic flow, where 
there is a non-zero vorticity which is slowly varying. The approach is to transform the 
equation of motion to a reference frame which rotates around the axis of vorticity. In this 
coordinate system, the strain field oscillates in directions which are perpendicular to the 
vorticity vector, in addition to having random temporal fluctuations. The effect of these 
oscillations is to reduce the effective intensity of the random strain field in directions 
perpendicular to the vorticity vector, so that an isotropic problem with vorticity is 
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transformed to an axisymmetric model with a velocity gradient which is a pure strain. 
This reduction was also discussed in [8], but is included here for the convenience of the 
reader. 

In order to isolate the effect of the vorticity in the equation of motion for the 
monodromy matrix, M, we introduce another monodromy matrix Mq which evolves 
under the vorticity alone: 

M = (S + n)M , Mo = Mo . (10) 

Note that Mo(t) is just a rotation matrix, describing rotation about an axis in the 
direction of the vorticity vector uj. The two monodromy matrices may be related by 
writing 

M{t) = Mo(t)K(t) (11) 

where K(t) is an evolution matrix which describes the effect of the shear. An elementary 
calculation shows that K has the equation of motion 

K = cr(t)K (12) 

where 

a- = Mo ^SMo (13) 

is obtained from S by applying a time- dependent rotation. Consider the form of the 
matrix cr. In the case where the vorticity vector is frozen, and equal to Qq, the matrix 
Mq is a rotation matrix: Mq = exp(f2ot). Without loss of generality we can consider 
the case where the vorticity is aligned with the 2;- axis, with magnitude Q = 2uj, where 
u is the rotational angular velocity, so that Mo is a rotation matrix of the form 



Mo = exp(riot) = K{ujt) 



If the elements of S are S, 



f cos ut — sin cot \ 
sin ut cos ut 

the elements of cr are 



c 
s 



-s \ 
c 

1) 



(14) 



/ c^Su + s^S22 + 2csS_ 



12 



Sii 



^ - S^)Si2 + Cs{S22-Sii) 

12 



s^Sii + 0^822 - 2csS, 



\ 



cSu - SS23 
CS23 + 56*13 
S33 



.(15) 



-S^)Sr2 + Cs{S22 

cSi3 - SS23 CS23 + 55*13 

Note that all of the off-diagonal components oscillate with angular frequency u or 2uj. 
The diagonal component in the direction of the vorticity vector does not oscillate, but 
the other diagonal elements contain both oscillatory terms and non-oscillatory terms. 



3.2. Limit of short correlation time for strain rate 

Now consider the case where the strain rate S is sufficiently small that the strain which 
accumulates over its correlation time Tg is very small. In this case the evolution of the 
matrix K (defined by equation (fT2l) ) can be described by a diffusive process. Specifically, 
we consider the evolution of ( fT2l) over a time period 6t which is large compared to the 
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correlation time of the strain field t^, but still sufficiently small that the strain which 
accumulates over this time interval is small. We write 

K{t + 6t) = (I + 6'E{6t, t)) K{t) (16) 

where the 51] are small and may be assumed to be random matrices, chosen 
independently at each timestep. We characterise the evolution ( |T2l) by computing the 
statistics of the random strain increments 6H, which are in turn obtained from the 
random strain S{t) using ( JT2l) and ( JT5l) . The advantage of considering the small elements 
S is that they are small random increments which are applied independently at each 
timestep. This enables their effect to be analysed using a Fokker-Planck equation. First 
consider the relation between the elements of the matrices ST, and cr. By integrating 
( !T2|) and using the definition ( !T6|) we obtain 

ft+St 

si:{6t,t) = j^ dt' (T{t') {i + 6i:{t' -t,t)) . (17) 

Iterating this expression, taking t = 0, and suppressing the initial time t in the arguments 
of 5 S we obtain 

(5S(5t) = / dti o-(ti) + / dti dt2 (T{ti)(T{t2) + 0{a^) . (18) 
Jo Jo Jo 

Using the fact that the correlation time is assumed to satisfy 6t ^ Ts, we can write 

6^{6t) = / dt cr(t) + ^ / dt (o-(t)cr(O)) + 0{6t^/^) . (19) 

Jo 2 J-oo 

The first of term is a random variable with mean zero and size 0((5t^/^), giving rise to 
a diffusion term in a Fokker-Planck equation. The second term represents a drift at 
a velocity which is well-defined in the limit as 6t — )■ 0. The remaining terms may be 
neglected. In order to formulate the Fokker-Planck equation, we must determine the 
statistics of the increments 5Sy(5t). 

If UTs ^ 1, the effect of the oscillatory terms in equation ( 1T5|) is negligible. Let us 
consider how to treat the problem when uTg is not small. To simplify the discussion, 
consider the quantity 

SF= dt /(t) cos(wt) (20) 
Jo 

where 6t/Ts ^ 1, and where /(t) is a random function which satisfies 

(/(t))=0, {f{t)f{t'))=C{t-t') . (21) 

The spectral intensity /(z/) of the fluctuations of /(t) is defined by 

I{u) = / dt exp(iz/t) C(t) (22) 



and we shall assume that C{—t) = C{t), so that /(— w) = I^uj)- The expectation value 
of SF is equal to zero. Its variance is 

fSt fSt 

{6F')= dtj dt2 (/(ti)/(t2))cos(a;ti)cos(c^t2) 
Jo Jo 

I rSt rSt 

= - dti / dt2 C{ti - t2)[cos(u;(ti - ts)) + cos(w(ti + ts))] 
/ Jo Jo 
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^6t[I{uj) + I{-uj)] + 0{6t^) = ^6tl{u) + 0{6t^) 



(23) 



The third steps assumes that u6t ^ 1, as well as 6t/T ^ 1. 

Now consider the effect of the random strain model defined by ([6])- ([9]) in the limit 
where the timescale Tg of the fluctuations of Sij(t) is very small. We assume that the 
functional form of the spectral intensity of each component Sij is the same, but that 
their variances are different, so that the spectral intensity of Sij{t) is {Sfj)!^^), implying 
that the intensity function is normalised so that /(O) = 1. We represent the effect of 
the randomly fluctuating strain field described by ( ITSi) by an effective strain field with 
diffusive fluctuations. Note that 5t is assumed to satisfy 6t/T ^ 1, despite being 'small'. 
By applying ( j23l) . variance of is 



5t 



St 



dti / dt2 



-(1 + cos2u;ti)S'ii(ti) + -(1 - cos2cjt2)5'22(ti) + sai2ujtiSi2{ti] 



X 



^(1 + cos 2a;)f:2) 5*1 1(^2) + ^(1 - cos 2wt2) 5*22(^2) + sin 2a;t25'i2(t2) 



/■oo X 1 

= 5t dr -[2 + cos2u;r] (5ii(r)5n(0)) + -[2 + cos2a;r] (S22(r)522(0)) 

J —00 o o 
+ i cos(2a;r) (5i2(r)5i2(0)) + ^[2 - cos(2u;r)] (5n(r)522(0)) + 0(5^^) 

= |[2 + I{2uj)]{Sl,) + |[2 + I{2uj)]{Sl,) + ^-^I{2uj){Sl,) 

+ ^l[2-I{2uj)]{SnS22) +0{5e) 

Using the same approach, the full set of non-zero covariances of JSjj is 

5t r 



(24) 



= = - (2 + I{2uj)){Si;) + (2 - I {2uj)) {8^,822) + 2/(2u;)(5i^2) 



((5Sii(5S22) 



4 

6t_ 

T 

6t r 



[2 - I{2uj)){8l,) + (2 + /(2a;))(5n522) - 21{2uj){8l^) 
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I{2u){8l^) - I{2uj){8n822) + 21 {2uj) {8^,) _ 

5t/(0) [{28l,) + 2{8n822) 

(S?3) = (SL) =StI{u;){8!,) . (25) 

Finally, we must consider the mean values of the increments 6T,ij{6t). As an example, 
consider the evaluation of From the second term in the right hand side of (JT9l) . 

we have 



(5Sn) 



St ^ 

dt Y.^a,,{t)a,m 
i=i 



2 j- 



6t roo 

2 

^2 2 



- / dt c'{8u{t)8nm + s'{8u{t)822m 

Z J— 00 



+ ic' - S^)(5i2(t)5i2(0)) + C(5i3(t)5i3(0)) 
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[(1 + I{2uj)){Si,) + (/ - /(2a;))(5n522) 



+ 2I{2u){Sl,) + 2I{u){Sl,)] . (26) 

Only the diagonal elements of 5S have a non-zero contribution to the mean at 0{5t): 
we define velocity coefficients yUj as follows 

5t 



[I + I{2uj)){Si^) + (/ - /(2a;)) (^11^22) 



2I{2uj){Sl,) + 2I{uj){Sl,) 



{ST.22) = fJ'2St = fliSt 

(5S33) = f^sst = ^ U{si,) + 4(^11^22) + ^ntu){si,) 



(27) 



3.3. Uniaxial random strain in three dimensions 



In sections 13.11 and 13.21 we showed how an isotropic model with frozen vorticity and 
rapidly ffuctuating strain can be represented by an axisymmetric model where the 
velocity gradient is a pure strain cr. In the limit where the strain which occurs over 
the correlation time Xg is small, the effect of this strain is represented by a product of 
matrices I + 5S, where the small increments 5S are independently distributed at each 
timestep of size 5t. They have diffusive ffuctuations, so that 55] = 0(5t^/^). The matrix 
cr is traceless, representing the fact that the velocity field is incompressible. The matrix 
(5X1 need not, however, satisfy tr(5S) = 0, although it is clear that the leading order term 
in (fT9l) is traceless. In this section we discuss how to parametrise such axisymmetric 
strain fields. 

We take this axis of rotational symmetry to be 63; the general case is obtained from 
this one by applying rotation matrices. The strain is described by a 3 x 3 matrix 
which we write in the form 

f 6A 5C SD \ f fiiSt \ 



5S 



6C 
6D 



V 








fii6t 

/i35t J 



(28) 



SB SE 
SE -{6A + 6B) J 

where SA, 6B, 6C, 6D and SE are random variables with mean value zero, and diffusive 
ffuctuations: {SA) = and {SA^) = 2DAASt, (SASB) = 2DABSt, etc. 

Applying a rotation about the axis by angle 6 to the random component of STi 
gives a transformed matrix, with elements SA', SB', SC , SD' and SE' , given by 

SA' = cos^ eSA + sin^ eSB + 2 cos sin BSC 
SB' = sin^ eSA + cos^ OSB - 2 cos ^ sin BSC 
5C' = (cos^ e - sin^ e)SC + cos 6 sme{SB - SA) 
SD' = coseSD + sin eSE 

SE' = cos eSE - sin OSD (29) 
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where c = cos^ and s = sin 6. The non-random diagonal component is invariant under 
rotation about 63. Note that 6A' + 6B' = 6A + 6B, so that the element is invariant 
under rotation. 

We require that the statistics of the elements are invariant under the rotation angle 
6. It is clear that 6 A and 6B must have the same variance, as must 5D and 5E. Without 
loss of generality, we can consider a model with {SA"^) = 26t. We therefore characterise 
the model by the following statistics, where a, (3, 7 are three constants: 

{6A^) = {6B^) = 26t 
{5A5B) = 2a6t 

{6C^) = 2f36t 

{6D^) = {6E^) = 2^6t . (30) 

Other covariances, such as {6B6E) , are equal to zero. The requirement that the statistics 
of the rotated matrix are independent of 6 leads to the equations 

{6A'^) =2[c^ + s^ + 2^s^a + Ach^l3]6t = 26t 
{A'B') = 2[-A^s^(3 + 2ch^ + (c^ + s^)a]6t = 2aSt 

(C"2) = 2[(c^ + s^- 2ch^)(3 + ch\2 - 2a)\6t = 2pSt . (31) 
Rotational invariance therefore leads to an equation which must the satisfied by a and 

a + 2(3 = I (32) 

so that the model for a uniaxial random strain has four independent parameters, which 
we can take to be a, 7, /xi and fi^. 

For a special choice of these parameters the model is isotropic. Clearly this 
requires fii = /i3, and (5C^) = {SD'^) = {SE"^), implying j = 13. Also, requiring 
{{SA + 5Bf) = {dA^) = {6B'^) gives 2 + 2a = 1. Solving these equations we find that 
the the covariances of the random terms are fixed in the isotropic case 

a = ~ , /3 = 7 = ^ , /i3 = /ii . (33) 

Another notable limit of the model is the case where the matrix is diagonal: this model 
is /3 = 7 = 0, implying a = 1. 

4. Alignment in random strain fields 

4.I. General solution in a diffusive axisymmetric strain 

In section [3] we described the construction of a model for the alignment of microscopic 
rods with vorticity, in which the velocity gradient is represented as a strain field with 
diffusive fluctuations, axisymmetric about the direction of the vorticity. First we 
consider the alignment of rod-like particles under a succession of independent random 
shears I + which satisfy the conditions derived in section [373] for the shear statistics 
to be uniaxial, before discussing the specific model for rod alignment in section 14.21 
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Using the approach summarised by equations (jl]) and (jSD, the direction vector n of 
a rod-hke particle evolves according to the linear equation 

(I + 5S)n(t) = (1 + 6R)n{t + 6t) (34) 

where ST, is the infinitesimal strain in time 6t, previously introduced in equation ( !T7|) . 
and 6R is the fractional change in length of the vector under the linear evolution 
equation. Write n{t + 6t) = n(t) + 6n + 0(5n^), where 5n ■ n = 0. Because of 
rotational symmetry about the 2;-axis, we can assume without loss of generality that 
the y component of n is equal to zero. We therefore consider the following orthogonal 
basis of unit vectors 

n = (sin6', 0, cos6') = {x,0,z) 
m = (— cos 9, 0, sin 9) = {—z, 0, x) 

k =(0,1,0) . (35) 
where 9 is the polar angle, and z = cos 9. Writing 6n = 6Xm + SYh, we have 

n(t + 6t) = n + 6Xm + 6Yk - ^{6X^ + 6Y^)n + 0{6n^) . (36) 

By taking the dot product of ( IMl) in turn with n, m and k, we find, respectively to 
leading order 

6R^n- 61211 = SEnn (37) 

and 

m-6n{l + 6R) ~ m-5Sn = 5Smn 

k-6n{l + 6R) ~ k-5Sn = 5Sfc„ (38) 

Let us characterise the evolution of n through the evolution of its projection onto the 
es axis, namely 

^ = es ■ n . (39) 

This is a convenient choice because z will have a uniform probability density function 
for an isotropic strain field. Using ( l36l) . we find that 

z + 6z = e3-n{t + 6t) = cos9 + sin 96X - ^ cos 9{6X^ + 6Y^) . (40) 
We define the drift velocity and diffusion coefficient of z by 

{6z) = vjt , {6z^) = 2DJt . (41) 
Using (142]) and dH?]), (IHHD we obtain 

vjt = - 5S„„5S^„) - ^{6Sl^ + 5SL) + Oi6t'/') (42) 



and 



DJt=^{l-z'){SElJ + 0{Sty'). (43) 
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Now consider that statistics of the fluctuations of z for the uniaxial strain model. For 



the model defined in section [3T3l we have 

5S„„ = dAx"^ + 26Dxz - {6 A + 5B)z^ + ^nx^St + /ig^^^t 

= 5D{x^ - z^) - {25A + 5B)xz + (^ug - ^i)xz5t 
5^kn =6Cx + 6Ez (44) 



where x = \/l — z^. We can combine these relations with (H2|) and (143!) to determine 
Dz and v^. 

I- z^ 

DJt = -^{[5D{1 - 2z^) - (25 A + 5B)xz]^) 

vjt = -x{[6A{l - 2z^) - 5Bz^ + 26Dxz][6D{l - 2z^) - {2b A + bB)xz\) 

- 2z^) - {2b A + bB)xzf) - ^{[5Cx + 5Ezf) + A/ix^z^t (45) 

where A/x = fis — fii- Using the statistics of the elements 6A, 6B, 6C, 6D and 6E, and 
ordering the resulting expressions as polynomials in z, we have: 

D, = ^{1- x^) [7 + (5 + 4a - 47)^2 _ (5 + 4« _ 4^)2^ 

7 5 5 , \ / 37 29 15 \ , 

4 + 4«--7 + A/.jz+(^----a + -7-A/.j 

/15 \ 

+ f y + 6a - 67] / . (46) 
The steady-state probability density for z, namely P{z), satisfies 

v.{z)P{z) = -^jDz{z)P{z)) . (47) 

In the isotropic case, we have a = —1/2 and 7 = 3/4. In this case we find 

3 3 
Dz = -(1 — z'^) , Vz = —-Z , (isotropic case) (48) 
8 4 

and the normalised solution is P{z) = |for— l<z<l. 

In the general case, we find that (1 — z^) is a factor of Vz — D'^, and the differential 

equation (HTll is 

1 dP _ [6(5 + 4a - 47)^2 _ 13 _ ^ ^07 + 4A/i] 

Pd7~ 4 [7 + (5 + 4a - 47)^2 - (5 + 4a - 4'y)z^] ^ 

it us useful to change the variable to u = z'^. In terms of u, the differential equation 
fH9|) may be written 

1 dP _ _6(5 + 4a - 47)7/ - 13 - 11a + IO7 + 4A/i 

P dw ~ 8 [7 + (5 + 4a - 47)^ - (5 + 4a - 47)^2] 
Representing the right-hand-side using partial fractions, we obtain 

1 dP _ c+ ^ c_ 

P d-u — M u — U- 

where u± are the roots of the denominator on the right-hand-side of ( l50l) 



(50) 
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and where the coefficients are 

_ (4A/i - 2a + 7 - 2)u± - 13 - 11a - 27 - 4Afi 

~ 4(5 + 4a) 

The probabihty density expressed in terms of z is then 

P{z) = C{z'^ - 

where C is a normahsation constant. 



-{z^-u^r^ 



(53) 



(54) 



4-2. Solution of rod alignment model 

Now we apply the solution obtained in section 14.11 to the model for alignment of 
microscopic rods, as developed in sections [2] and [31 In section 12.21 we introduced the 
Ornstein-Uhlenbeck model for a random, isotropic velocity gradient field. The theory in 
section [3] made two assumptions. In section [3TT] it was assumed that the vorticity varies 
slowly, and section 13.21 made a further assumption that the strain field is small. Let us 
consider the implications of these assumptions for the parameters of the model. The 
assumption that the vorticity varies slowly implies that is large compared to other 
timescales in the system of equations. The typical strain rate IS"! = 



'(tr(S2)) and the 

correlation time Tg should satisfy \S\ts <^ 1. The solution of the Ornstein-Uhlenbeck 
process imphes 

(tr(S2)) = IOAts (55) 

so that the criterion for the strain to be small is simply D^t^ ^ 1. The angular velocity 
uj is related to the magnitude of the vorticity f2 by i7 = 2uj. The magnitude of the 
vorticity is estimated by (p?) = — |tr((f2^)) = SD^Tv The rotation rate u has a 
Gaussian distribution, with variance 



cr 



(56) 



Because the Ornstein-Uhlenbeck model has an exponential decay of correlations, given 
by equation (|8]), the spectral intensity of the strain fluctuations is a Lorentzian function: 

1 



(57) 



1 + z/2r2 

In order to apply the results in section 14.11 we must specify the covariances of the 
fluctuations of the axisymmetric effective strain tensor. If, in accord with the notation 
of section [3731 we normahse the variances so that {Sf^) = 1, (5'ii5'22) = a, {8^2) = /3, 
(S'^g) = 7, the non-zero covariances and expectation values of are 

{5El) = {SEl2)=6t^-- ^ 
((5Sii(5S22) 
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(S?3) = (E^3) = 6tl{ujh 
{6El,) = 2 St I {0){1 + a) 

((5Sn) = {6^22) = St 
(5S33) =6t[I{0){l + a) + I{u)-f] . (5^ 



/(0)(^ + f)+/M| + /(2u; 




We use the assumption that the original random strain field Sij is isotropic, so that the 
statistics of these elements satisfy fl33l) . Using fl58|) we obtain 

(5S?,) = (5Sy = ^[l + 3/(2a;)] 
(5Sn5S22) =^[l-3/(2u;)] 
^2) = j3/(2a;) 

(5E?3) = (SEl,) = ^3/(u;) 

(5S^3) = 6t 

((5Sn) = = f [1 + 3/(c^) + 6/(2a;)] 

(5S33) = ^ [2 + 3/(0;)] . (59) 

Normalising these by dividing by the anisotropy of the strain field induced by 

the vorticity is characterised by modified forms for the parameters defined in (130|) and 
a scaled value of the parameter A/i appearing in (H6l) : 

1 - 3/(2a;) 



a 



1' 



1 + 3/(2u;) 

3/(2a;) 
1 + 3/(2w) 

3/(0;) 
1 + 3/(2w) 



_ 3H-/M- 2/(2.) 

^ 2 l + 3/(2u;) ^ ^ 

Note that a' + 2/3' = 1, as expected from fl32l) . Because J(0) = 1, when the vorticity is 
zero, we have a' = — 1/2, = 7' = 3/4, so that we recover the statistics of an isotropic 
strain field. In the limit as w —t- 00 we expect I{2uj) <^ I{uj) <^ 1, so that to leading 
order we may set a = 1, /3 = 0, A/i = | and regard 7 as a small parameter. 

Let us consider how to evaluate the probability density function of 2; = n ■ in 
the limit where lot^ ^ 1. The general expression for the probability density is equation 
flM|) . In this limiting case, we may approximate the coefficients a and 7 by 

a ~ 1 , 7 ~ ' ~ ^ • ^^^^ 
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Figure 1. Comparison between PDF of z = n • 63 obtained by simulation of Jeffery's 
equation of motion for the random strain model, and the theoretical prediction, 
equations (f52|) - (|54[) and ([60|) . In these simulations the vorticity is frozen so that 
C = wTs = 1,3,5,7. 



When 7^1, the poles u± of the probability density function and the coefficients c± in 
( l54l) are approximated by 



1 

1 
2 



1 



1 
9 



c 1 ~ 



-1 



so that the probability density function is approximated by 

-1/2 / rs, \ -1 



Pjz) -C{z' + 1 



1 + ^ 
9 



(62) 



(63) 



where C is a normalisation constant, and where the subscript a; is a reminder that this 
distribution is evaluated for a fixed value of u. 

We verified this relation by simulating the orientation of rod-like particles using 
equations and ([5]), using the Ornstein-Uhlenbeck model for the velocity gradients. 
The components of the vorticity were frozen, so that the only non-zero elements are 
fli2 = —fl2i = The PDF of 2; = n ■ 63 is plotted in figure [T]for = uTs = 1, 3, 5, 7, 
showing good agreement with the theoretical prediction, equations ( l60l) and ( I52|) -(l5^. 
The numerical simulations used = 1, Ds = 10~^, and the timestep of the numerical 
integration was dt = 10 ~^ or smaller. 
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Figure 2. Theoretical PDF of z = n • 63 for the model where the vorticity is 
frozen, with magnitude = 2a;. The distribution is plotted for three values of 
C = wTs, namely ^ = 1,3, 10. By comparison we have also plotted the exact probability 
density obtained from equations ([52 l) -(f54 l) and (|60| . At C = 10 the different plots are 
indistinguishable. 



In figure [2] we plot the theoretical PDF of z for three different values of the 
dimensionless variable ( = uTs, comparing (!63l) with the exact expression obtained from 
using fl60|) in equations f l52|) - fl5^ . For uTg = 10 plots of the exact and approximate 
PDF lie on top of each other. We observe that as ur^ — ?■ 00 the distribution becomes 
concentrated around z = ±1 (rods perfectly aligned with the vorticity) and around z = 
(rods ahgned perfectly perpendicular to the vorticity vector). The peak at 2; = ±1 is 
seen to be higher but narrower. In figure[3]we plot {\z\) and {z'^)- Both of these statistics 
approach | in the limit as C — ^ cxd, indicating that in this limit both peaks carry half of 
the probability. 

In practice the magnitude of the vorticity, u, is not frozen but fluctuates slowly. It 
has a Gaussian distribution, with a variance cr^ = (w^) = ^D^t^. Our final estimate for 
the probability density of z is, therefore, the result on integrating the normalised PDF 
given by f l63|) over a Gaussian distribution of u: 

2 

P{z) = / da; expi-tuy2a^)PUz) . (64) 

V 27rcr Jo 

This PDF depends upon a single dimensionless parameter 1^ = arg. The functions 
obtained by numerical evaluation of the integral in f lM|) are plotted in figure H] for three 
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Figure 3. Moments (|z|) and (z^) as a function of = ujts. 

different values of ,^ = aTg. We used the exact formulae for Puj{z), equations (l52l) - (!54ll 
and fl60|) . because the integral includes the region where uTs is small. 

5. Discussion 

We have determined the distribution of z = n-e^^ analytically for a model of microscopic 
rods in a random velocity field with isotropic statistics. The PDF shows a maximum 
at z = 1 corresponding to alignment parallel to the vorticity, similar to findings of DNS 
studies of Navier-Stokes turbulence |8]. 

We conclude by making a few remarks about the relationship between the regime 
which we have analysed and the velocity gradient statistics for Navier-Stokes turbulence. 
The model for the velocity gradient of an isotropic flow which was introduced in section 
12.21 has four parameters, namely t^, Ts, and D^, all of which have dimensions which 
depend only upon time. There are, therefore, three dimensionless parameters. In our 
analysis the vorticity was frozen, so that Tv — )■ oo. The magnitude of the vorticity, which 
is of order u ~ a/IVtv was chosen so that C, = ut^ is finite. The diffusion coefficient 
was assumed to be very small, so that the fluctuations of the strain are very small and 
may be treated using a Fokker-Planck equation. 

In fact the form of the Navier-Stokes equation restricts the choice of parameters 
in the Ornstein-Uhlenbeck model for the velocity gradient: it is well known that 
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Figure 4. Theoretical PDF of z = n • 63, averaging over tire slowly varying 
vorticity parameter ui, which is Gaussian distributed with variance ct^ = jD^t^. The 
distribution ([61]) is plotted for four values of ^ = ctts, namely ^ = 1.15, 3, 5, 7. 



tr(02) +tr(S2) = [16J, which gives a further relation between Ds and D^. The Navier- 
Stokes equation also implies that the rate of dissipation per unit mass is £^ = i/tr(A'^A), 
which enables the norm of the velocity gradient to be expressed in terms of the 
Kolmogorov time, tk = 



'u/S. These results imply the following relations, which 
determine the ratio of the diffusion coefficients Dg and (see [8]): 

^ " ^ (65) 



Numerical studies indicate that the exponential correlation function is a reasonable 
model for the statistics of fully developed turbulence, with the parameters Tg and 
satisfying Tg ~ 2.3tk and ^ 7.2tk (these are the values quoted in |8], which discusses 
earlier work on the velocity gradient correlation functions). This justifies the assumption 
that the vorticity is slowly varying, and the variance of the vorticity parameter is 
estimated to be cr^ = (cj^) = f-Dv^v, so cr = 1/{2tk), implying that ^ = aTs ~ 1.15 
is the value which should be compared with the data on alignment in Navier Stokes 
turbulence, discussed in [8]. There is a qualitative but not quantitative agreement 
between the curve in figure H] for ^ = 1.15 and the results of DNS simulations in figure 
2 of [8j: both show a peak in the PDF at z = 1, but this peak is more pronounced 
in the DNS data. We conclude that our model should be understood as a laboratory 
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for understanding alignment of microscopic rods with vorticity, rather than providing a 
quantitative description. 
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